Simple and Robust Solver for the Poisson— Boltzmann Equation 
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A variational approach is used to develop a robust numerical procedure for solving the nonlinear 
Poisson-Boltzmann equation. Following Maggs et al., we construct an appropriate constrained free 
energy functional, such that its Euler-Lagrange equations are equivalent to the Poisson-Boltzmann 
equation. This is a formulation that searches for a true minimum in function space, in contrast to 
previous variational approaches that rather searched for a saddle point. We then develop, implement, 
and test an algorithm for its numerical minimization, which is quite simple and unconditionally 
stable. The analytic solution for planar geometry is used for validation. Some results are presented 
for a charged colloidal sphere surrounded by counterions, and optimizations based upon Fast Fourier 
Transforms and hierarchical pre-conditioning are briefly discussed. 

PACS numbers: 02.60.Lj, 02.70.Bf, 41.20.Cv, 47.57.-s, 87.10.Ed 



I. INTRODUCTION 

In many soft-matter and biological systems electro- 
static interactions have a strong influence on the phys- 
ical behavior. A typical example are charge-stabilized 
colloidal dispersions [l| , whose structure is mainly deter- 
mined by the interplay between van der Waals attraction 
and electrostatic repulsion. A simplifying feature is the 
fact that it is often sufficient to describe the structure 
of the "cloud" of surrounding ions (counterions and salt 
ions) just in terms of Mean Field theory, in particular 
if the ions are monovalent. At the center of this theory 
is the well-known Poisson-Boltzmann equation, which, 
from the mathematical point of view, is a nonlinear par- 
tial differential equation. Analytical work in this field has 
been mainly confined to the linearized (Debye-Hiickel) 
version, which however is often insufficient, in particular 
near strongly charged objects. This has prompted efforts 
to develop methods to solve the equation numerically, 
ideally in three dimensions and without restrictions on 
the underlying geometry or spatial symmetry. These are 
usually based on standard finite-difference Q or finite 
element Q techniques (for a recent review see Ref. Q), 
and have meanwhile reached a substantial degree of so- 
phistication and complexity. 

Recently, however, Maggs Q has put forward a com- 
pletely new approach to electrostatics in soft -matter re- 
search, which bears a certain similarity to lattice gauge 
theories @. The central idea is to use the electric field 
E instead of the electrostatic potential t/j as the quan- 
tity on which the algorithm operates, and to view Gauss' 
law eV • E = p as a constraint for the field configura- 
tions. On a lattice, it is then easy to construct an elec- 
tric field (E) configuration which satisfies Gauss' law for 
a given (arbitrary) charge density p. The hard part is 
rather the transversal part of the field, which should sat- 
isfy V x E = 0, but initially does not (unless one uses a 
sophisticated initialization procedure based upon solving 
the Poisson equation). This transversal degree of free- 
dom can then be removed by local relaxations (this is 



the approach taken in the present paper), or integrated 
out by performing a Monte Carlo p, 0] or Molecular Dy- 
namics [1, [t| [l(| simulation on the overall system. A 
crucial aspect of the method is to locally update both p 
and E simultaneously in a way that Gauss' law is still 
satisfied after the update, i. e. the "constraint surface" 
is never left. 

The original Maggs approach is based upon a system of 
discrete charges whose statistical physics is treated in a 
consistent way. However, it is also possible to apply this 
to the Mean Field version of the theory, i. e. the Poisson- 
Boltzmann equation. The purpose of the present paper 
is to outline how this can be done in practice. In Sec. 
HI1 we derive the method by re-formulating the Poisson- 
Boltzmann theory in terms of a free energy functional, 
which is minimized by using a Maggs-type algorithm. 
In contrast to previous formulations, this functional pro- 
vides a true minimum, and therefore the procedure is 
very straightforward and simple, unconditionally stable, 
and has a rather modest storage requirement which scales 
only linearly with the number of grid points. In terms 
of computational speed, the method can probably not 
yet compete with the existing packages; however, it is 
reasonable to assume that more advanced versions that 
combine the basic methodology with acceleration tech- 
niques like adaptive mesh refinement and / or unstruc- 
tured meshes, may become a very useful tool. Section 
IIIII presents numerical results obtained with our current 
simple implementation, while Sec. II VI discusses optimiza- 
tions based upon Fast Fourier Transforms and hierarchi- 
cal pre-conditioners. Finally, Sec. fVl finishes with some 
concluding remarks. 



II. DERIVATION OF THE ALGORITHM 

A. Poisson-Boltzmann Equation 

Consider a system of fixed charges, with total charge 
Ze (e > denotes the elementary charge) dispersed in a 
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solvent with dielectric constant e. The system is confined 
to a three-dimensional finite domain of volume V. The 
charges can be distributed either in macroscopic particles 
or in any other kind of boundary. Furthermore, the do- 
main contains counterions of total charge —Ze such that 
the system as a whole is charge neutral. The total num- 
ber of counterions is No, such that Z = —zqNq, where zq 
is their valence. Furthermore, the presence of other ionic 
species (salt ions) is allowed if they satisfy charge neu- 
trality, i. e. if J2a>i z aN a — 0, where z a is the valence 
of the ionic species a and N a is its number. 

The continuum theory provides equations for the con- 
centrations (particle number densities) c a (r) of each ionic 
species a and the electrostatic potential ip(r). In the 
stationary regime (no time dependence) and in thermal 
equilibrium characterized by the thermal energy ksT, 
they take the form [l[ 



fcgTV lnc Q + ez a Vip = 0, 

sV 2 lp + 2_. ez aC a = 0. 



(1) 

(2) 



The first equation is the equilibrium version of the 
Nernst-Planck equation, which balances the diffusion 
current of ionic species a against the drift caused by the 
electric field E = —Wip. The second equation is the Pois- 
son equation, taking into account all ionic species as a 
source term, while the fixed charges appear as boundary 
conditions. The total number of particles is obtained by 
integrating the concentration over the whole domain: 



c a dV, 



(3) 



and the charge neutrality may then be expressed by 



N j z a I c a dV = —Z, 
Jv 



(4) 



which is an equation for ip only. However, for the algo- 
rithm to be discussed below, it will be advantageous to 
rather consider the equivalent coupled set of equations 
(Eqs.HlO. 



B. Reduced Units 

The Poisson-Boltzmann equations can be rewritten in 
terms of nondimensional quantities, i. e. in reduced units. 
The argument of the exponential in Eq. [5] suggests the 
most natural way of rescaling the potential: 



eip 



(9) 



i. e. the reduced potential is the electrostatic energy of 
an elementary charge in units of the thermal energy. This 
leads to 



Vlnc Q + z a Vip' = 0, 



(10) 
(11) 



where Ib = e 2 / {ATreksT) is the Bjerrum length. In- 
troducing a parameter characteristic length 
scale (see below), the gradient operator is rescaled via 
V = kV'. This allows writing the equations in the nondi- 
mensional form 



where 



W + E^ = o, 



(12) 
(13) 



(14) 



with 



y^Zq / CadV = 0. 
JV 



By integrating Eq. [T] one obtains 

ez a tp 



A a exp 



(5) 



(6) 



where the integration constant A a has the dimension of a 
concentration, and, for normalization reasons, must have 
the value 



N„ 



J v exp (—ez a ip/kBT) dV 



(J) 



Inserting this result into Eq. [21 one obtains the Poisson- 
Boltzmann equation: 



eV 2 f/' + J]] cz a A a exp 



ez a ip 
~k^f 



0. 



(8) 



is the reduced concentration. In terms of the electric field 
E' = —V'tp', the equations are 



V In 4, 
V' • E' 



z a E', 



V'xE' = 0. 



(15) 
(16) 

(17) 



The normalization condition for the amount of species a 
is transformed to 



with 



c'dV = N' 



N' = A<Kl B nN a . 



(18) 



(19) 



The choice of the parameter k is completely immaterial 
for the mathematical formulation of the problem. It is 
only important to map the numerical results back onto 
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a physical system, and therefore a matter of convention. 
For many applications, the choice 



while variation with respect to E yields 



4ttZ, 



V 



(20) 



turns out to be quite useful: This is the natural screen- 
ing parameter in the finite-volume version of linearized 
Poisson-Boltzmann theory (therefore it also appears in 
the corresponding treatment of asymmetric electrolytes 
[TTI| ) . In the case of only one monovalent ionic species 
(the counterions a = 0), this reduces to 



V. 



k 2 =4tt2 b ^1. (21) 

It should be noted that in this case Nq 

From now on, we will be concerned with the problem 
of numerically solving the reduced set Eqs. [T5HT71 In 
what follows, the primes will be omitted, with the un- 
derstanding that all quantities (including N a and V") are 
given in reduced units. 



C. Variational Approach 

Following the ideas of Maggs and Rosetto @, the 
Poisson-Boltzmann equation can be re-formulated as 
a constrained variational problem, where a free energy 
functional is minimized. This functional is constructed 
such that its Euler Lagrange equations are equivalent to 
Eqs. EMU Its form is 



fdV, 



(22) 



/ = ^E 2 +^c Q lnc Q -i/> IV-E-^z Q c Q J 

a. \ a / 

- ( c «-ir) ■ 



(23) 



The first term corresponds to the electrostatic energy and 
the second to the entropy. Gauss' law and the mass nor- 
malization conditions are included as constraints, via La- 
grange multipliers: The field ip (r) is the electrostatic po- 
tential, while the numbers /i Q are the chemical potentials 
of the species a. 

The Euler-Lagrange equations for this variational 
problem will be derived next. Variation with respect to 
■0 and fi a just recovers the constraint equations 



V • E = y^ZqCq, 

a 

c a dV = N a . 

Variation with respect to c a results in 

lnc Q + 1 + Tpz a — fj, a = 0, 



(24) 
(25) 

(26) 



E 



-VV>. 



(27) 



The dependence on the unknown Lagrange multipliers is 
removed by taking the gradient of Eq. [521 and the curl of 
Eq. [27] to obtain 



V lnc a + z a Vip = 0, 
V x E = 0. 



Finally, inserting Eq. [27] into Eq. [35] leads to 
Vlnc Q = z a E. 



(28) 
(29) 



(30) 



In summary, the derived equations (Eqs. [30] (2H [55]) are 
the desired set — in other words, the solution of the 
minimization problem is identical to the solution of the 
Poisson-Boltzmann equation. 

In this context, one should make a few important 
observations. Firstly, we note that solving the Euler- 
Lagrange equations will provide one and only one solu- 
tion: It has been proven (by variational methods) that 
the Poisson-Boltzmann equation has one unique solution 
[l2T [l3| . Secondly, it is easy to show that the solution is 
indeed a local minimum: By writing 



E 



E<°> + SE, 
c^+6c a , 



(31) 
(32) 



where E^ ^ , c a form the solution of the problem, and SE, 
Sc a are small deviations which satisfy the constraints, i. e. 



V-SE = ^2z a Sc a , 

a 

Sc a dV = 0, 



(33) 
(34) 



one finds 



F = 



iE(°> 2 + ^ci )ln C i°)j dV 

a c a ) 



i. e. small deviations from the solution will always in- 
crease the functional. Together with the uniqueness, this 
shows that the functional has one and only one mini- 
mum, which corresponds to the solution of the Poisson- 
Boltzmann problem. Therefore, a procedure that relaxes 
all degrees of freedom of the functional such that it is 
systematically decreased, while staying on the constraint 
surface, will ultimately run into the one and only min- 
imum of the free energy landscape. The iterative pro- 
cedure may be slow and hampered by small eigenvalues 
of the Hessian at the minimum, but such problems can 
be kept under control by careful convergence checks and 
variation of the number of iterations. A simple algorithm 
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that initializes the system on the constraint surface, keeps 
it there, and systematically decreases F by local updates 
of all non-constrained degrees of freedom will be outlined 
in Sec. Ill El Since the constraints are always satisfied 
during the procedure, the Lagrange multiplier terms in 
the functional may be omitted, such that it is simplified 
to 



F = 



^ j^E 2 + ^c Q In c Q j dV. 



(36) 



It should also be noted that previous approaches that 
were also based upon a free-energy functional (see, e. g., 
Ref. [Hj]) , did not truly search for a minimum, but rather 
for a saddle point. These methods were not based upon a 
functional that involves the electric field, but rather one 
that depends on the electrostatic potential, 



F 



fdV 



sr i - HA 

/ , \ c a y J , 



(37) 



(38) 



where the before, the Lagrange multipliers cor- 

responding to the mass normalization conditions. How- 
ever, in contrast to its meaning in Eqs. [22l [23} ip is 
here not a Lagrange multiplier, but rather a degree of 
freedom. It is straightforward to show that the Euler- 
Lagrange equations of this problem are equivalent to the 
Poisson-Boltzmann equation. To show that this solution 
is indeed not a minimum but rather a saddle point, one 
again decomposes and c a into the solution plus small 
deviations, 



- r (o) 



V» (0) + Sip, 
5c a , 



with 

This yields 
F 
fo 



Sc a dV = 0. 



(fo + h) dV + O (6c 3 

r 

-^( V ^ (0) ) 2 + E c i 0)ln ^ 

a 

^ (0) E^ 0) < 

a 

a c a 



(0) 



(39) 
(40) 



(41) 

(42) 
(43) 

(44) 



i. e. the quadratic form of the deviations is not positive- 
definite. It is natural to suspect that this lack of positive- 
definiteness causes various numerical difficulties in terms 
of stability, which therefore are intrinsically absent in the 
new formulation. 



D. Discretization 

The computational domain is a rectangular paral- 
lelepiped of size l\ x I2 x Z3 with periodic boundary con- 
ditions. This box is discretized by a simple orthorombic 
(usually: cubic) lattice with sites r and lattice spacings 
Axi, i = 1,2,3 enumerating the Cartesian directions. 
The volume of a unit cell is thus AV = Ax± Ax2 AX3. 
The concentrations c a are variables on the sites, while the 
electric field is associated with the links. The positions 
of the concentration fields are the vectors 



r (n) = (Axinx, Ax 2 n 2 , Ax 3 n 3 ) , 



(45) 



where Uj are integers. The field E\ is located at the 
positions 



ri(n) = (Axi(ni + 1/2), Ax 2 n 2 , Ax 3 n 3 ) . 
Similarly, the positions for E 2 and E 3 are 

r 2 (n) = (Ax 1 n 1 ,Ax 2 (n 2 + l/2),Ax 3 n 3 ) 



(46) 

(47) 
(48) 



r 3(n) = (Aaiini, Ax 2 n 2 , Ax 3 (n 3 + 1/2)) . 

respectively. Furthermore, it is useful to define 

ri(n) = (Ax 1 (n 1 -l/2),Ax 2 n 2 ,Ax 3 n 3 ) , (49) 
r 2 (n) = (Ax inu Ax 2 (n 2 - l/2),Ax 3 n 3 ), (50) 
r 3 (n) = (Ax 1 n 1 ,Ax 2 n 2 ,Ax 3 (n 3 -l/2)). (51) 

These definitions allow to approximate the functional by 

C2) 



F 



n i=l 

+ ^c a (r (n))hc o (r (n)) : 



and to also discretize the divergence operator in a 
straightforward way: 

(V • E) (r (n)) = £ -L (25< (r,(n)) - E t (r{(n))) . 



(53) 



Gauss' law then reads 
3 



£^(^(11))-^ (1^(11))) 



Ax. 

Introducing fluxes via 

01 = 

02 : 

03 : 



£z Q c a (r (n)) 



EiAx 2 Ax 3 , 
E 2 Ax 3 Ax ly 
E 3 Ax 1 Ax 2 , 



(54) 

(55) 
(56) 
(57) 
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this is rewritten as 

3 

]T (n(n)) - & (r^(n))) = AV £ (r (n)) . 

i=l a 

(58) 

Finally, the normalization condition for the amount of 
ionic species a is discrctizcd as 

5>(r (n)) = |^. (59) 

n 

E. Algorithm 

The numerical minimization procedure starts from 
some configuration of the discretized fields c a and E 
which satisfies all the constraints, i. e. the normaliza- 
tion conditions for the ions, plus Gauss' law. The al- 
gorithm then performs successive local changes in the 
electric fields and the concentrations, analogously to the 
Monte Carlo moves of Maggs and Rosetto ||. These 
moves have the big advantage that they rigorously con- 
serve the constraints. In contrast to Ref. [5J, however, 
the moves are not stochastic, but rather deterministic, 
and constructed in such a way that they decrease the 
functional, and do this optimally. Since the free energy 
landscape of this problem has a simple structure (as dis- 
cussed in Sec. Ill C[) . the procedure relaxes the fields into 
the one and only minimum, which is the solution of the 
Poisson-Boltzmann equation. 

The algorithm can be summarized as follows: 

1. Distribute the fixed charges. 

2. Classify the grid points. 

3. Distribute the ionic species uniformly in the move- 
able nodes. 

4. Initialize the electric field. 

5. Perform the field moves for all the plaquettes 
(smallest closed loops) in the grid. 

6. Perform the concentration moves for all pairs of 
adjacent moveable nodes. 

7. Check if the changes in the functional caused by 
steps 5 and 6 are less than a given tolerance: if yes, 
then stop, otherwise, return to step 5. 

In the beginning, fixed charges and ionic species must 
be distributed over the grid. Fixed charges are usually 
associated with surfaces. Therefore elements of surface 
charge density must be mapped onto elements of volume 
charge density, so that they can be associated with some 
of the nodes. These nodes are then marked as "fixed" 
and no particles can enter or leave them after initializa- 
tion. Further nodes may be marked as "fixed" if they 
are known to be empty (for example, if they represent 



the interior of a particle). The nodes representing the 
volume where ions can move are marked as "moveable" . 
Initially, each ionic species is uniformly distributed over 
the "moveable" nodes. Choosing the correct amount of 
charges then automatically results in charge neutrality of 
the overall system. 

The next step consists in initializing the electric field so 
that it satisfies Gauss' law for the initial charge distribu- 
tion. One possibility is to solve the Poisson equation with 
some numerical method. This needs to be done only once, 
in the initialization. An alternative, based on charge neu- 
trality, is to initialize each component by means of a re- 
cursion over the spatial dimensions [lol ]. which is equiv- 
alent to applying Gauss' law to linear chains of nodes. 
First the lattice is decomposed into a set of planes per- 
pendicular to the xi-axis and it is required that E\ takes 
the same identical value for all links with identical x\- 
coordinate. Then in one (arbitrary) plane of links we set 
E\ = 0. Starting from there, we can then calculate E\ 
step by step in the subsequent planes of links, where the 
change in E\ is given by the plane-averaged charge den- 
sity between the links. Assuming that the procedure is 
started at the charge plane X\ = 0, it reads 

E 1 (-0.5Ax 1 ,n 2 Ax 2 ,n 3 Ax 3 ) = 0, (60) 
Ei((ni + 0.5)Axi,n 2 Ax 2 ,n 3 Ax 3 ) = 
Ei((m - 0.5)Axi,n 2 Ax 2 ,n 3 Ax 3 ) + 
A Xl (p) (mAari), (61) 

where (p) is the plane-averaged charge density at x\ — 
n\Ax\. Charge neutrality combined with the periodic 
boundary conditions ensures that this procedure will give 
consistent results after closing the one-dimensional loop. 
Then each plane is decomposed into a sequence of lines, 
perpendicular to the x\- and the X2-axis, and the anal- 
ogous procedure is applied to obtain the field in x 2 - 
direction. The charges which occur here are the line 
averages, where however the plane averages have been 
subtracted (the latter have already been taken into ac- 
count via Ei). Finally, the lines are decomposed into 
sites, and £"3 is determined from the remaining charges 
where both line and plane averages have been subtracted. 

For field changes, elementary closed loops on the faces 
of the unit cells (plaquettes) are considered. For the or- 
thorombic lattice, these are comprised of four nodes and 
respective links, such that each node is connected to two 
plaquette links. Now, these four fields are modified in 
such a way that the flux on each link is changed by the 
same amount (taking into account the orientation along 
the closed loop). Therefore, Gauss' law will still be sat- 
isfied after that move, since at every node there will be 
some more flux entering but also the same amount of flux 
leaving. Let us, for example, consider a plaquette per- 
pendicular to the X3-axis, with a sequence of fields E\, 
E' 2 , E[, E 2 along the loop, where Ei, E[ are positive if 
the field points in positive xi-direction (and analogous 
for E 2l E'2). Then the field updates are given by 

Ei -» Et+SEk, (62) 
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E' 2 
E 2 

or, in terms of fluxes, 



E' 2 



SE' 2 , 



E[+SE[, 
Ei + 5E2, 



6<fii = Ax 2 Ax 3 SE\ = 5(f), 

6(f>' 2 = AxiAx 3 SE' 2 = S(f>, 

ctyi = Ax 2 Ax 3 SE[ = -84 

S(f> 2 = Ax 1 Ax 3 SE 2 = -64 



(63) 
(64) 
(65) 

(66) 
(67) 
(68) 
(69) 



where the parameter 60 can be chosen arbitrarily with- 
out violating Gauss' law. The associated change in the 
functional is given by 



SFAV 



{(A Xl ) 2 + (Ax 2 f) (S<f>) 2 

+ AVAxi (E 1 -E[)5(j) 
+ AVAx 2 (E 2 -E' 2 )54p, 

which is minimized for 

6<t> ~ 2 (A Xl f + {Ax 2 f 

x {Ax 1 (E[-E 1 )-Ax 2 (E 2 -E 2 )}. 



(70) 



(71) 



This yields the optimal values for the field changes. 

For concentration moves between two adjacent nodes 
(connected by a single link), Gauss' law is also conserved 
if the electric flux is updated accordingly. Suppose that 
before the move the two adjacent nodes and 
have, respectively, concentrations c^' and c^ of some 
ionic species with valence z. Without loss of generality 
we may assume that node B has a larger index value 
than node A. The electric flux from node A to node B is 
then given by (AV/ Al)E, where Al is the length of the 
link, and E the field on it. Now, a certain (positive or 
negative) amount Sc is moved from A to B, i. e. 



JB) 



JA) 
JB) 



Sc, 
Sc. 



(72) 
(73) 



Gauss' law tells us that the flux should change by 
— AVzSc, and this is the case for 



SE = -AlzSc. 
The resulting change in the functional is given by 

SE 



AF 
AV 



1 



E+^SE 



(74) 



(75) 



+ c^ In ( 1 



— 5c In 



6c 

c^ - Sc 



+ c^ In 1 



Sc 

7b) 



d B ) +8c 

In order to find the optimal value for 5c, we minimize 
this expression. The result is a nonlinear equation, 

C (A) _ C (B) cxp (_ Z AZ(£ - zAISc)) 



5c = 



1 + exp (-zAl(E - zAISc)) 



(76) 
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FIG. 1: Profile of the counterion concentration along the x- 
axis, for various grid resolutions as indicated by the legends. 



which must be solved numerically. By introducing 



^> + c(*>) 



JB) 



£ = -zAl (zAISc - E) , 



the equation is transformed to 

2 E 



tanh £ + 



(zAZ) c+ zAlc+ 



C- 

c+ 



0, 



(77) 
(78) 
(79) 

(80) 



which shows that it has exactly one solution (the slope 
of the left hand side is always positive). Since — 1 < 
tanh£ < 1, the solution will satisfy the condition 

2 E c 

-K 2 — £ t-. 1 < 1, (81) 

{zAlfc + zAlc+ c+ 

which is equivalent to c^ A ' — Sc > 0, c^ s ' + Sc > 0, such 
that Sc will be in the physically admissible range. Finally, 
the shape of the left hand side guarantees that a Newton 
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FIG. 2: Relative errors of the concentration profiles, for dif- 
ferent grid resolutions as indicated. 



iteration starting at £ = will always converge; hence, 
this rapid procedure was implemented. 



III. NUMERICAL RESULTS 

In this section, the feasibility of the present approach 
is demonstrated by two numerical examples. The choice 
of parameters is inspired by previous computer simula- 
tions on electrokinetics done in our group (lEI . [13]. 
We therefore quote them here in unsealed "physical" 
units, where Ao denotes our elementary length scale (the 
Lennard- Jones diameter in our simulations). All calcu- 
lations are done with a Bjerrum length lg = 1.3A . The 
fixed charge distribution (boundary condition) consists of 
positive charges only, while there is only one ionic species, 
the monovalent counterions with valence z = — 1. The 
resulting data are also given in "physical" units. 



A. Double Plane with Counterions 

Consider that the fixed charges are distributed in 
two infinite parallel plates, perpendicular to the x-axis, 
placed at x = — a and x = a. The surface charge den- 
sity in each plate is a. Furthermore, suppose that there 
are no salt ions and that the counterions (of valence z, 
concentration c(r)) are distributed in the region between 
the planes. Charge neutrality is given by 



/ ads+ i 
J s Jv 



zc(r)dV = 0. 



(82) 



For this case, the problem is one-dimensional and its 
analytical solution is well-known [l8j]. Therefore, it is 
ideally suited to test the algorithm. In reduced units, 
the Poisson-Boltzmann equation is 



d 2 ip 
dx 2 



-Az exp (—zip) 



(83) 
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Code B: 32 x 32 x 32 



O 



n 



Code A: 64 x 64 x 64 



Code B: 64 x 64 x 64 



o 



o 



„ 



Code A: 128 x 128 x 128 Code B: 128 x 128 x 128 
FIG. 3: Counterion concentration for a charged colloid. 



= —^r(-Aexp(-zip)), 

which can be interpreted as Newton's equation of motion 
of a particle with unit mass, whose coordinate is tp and 
where the time corresponds to x. Hence, this equation 
can be solved via standard methods of classical mechanics 
[lj|. The result is given by 



c(x) 



ip = — In cos (s— \ 

z \ a/ 

2 s 2 



where the parameter s is related to A via 2s 2 
The surface charge density is given by 



dx 



2 s 
tan s; 

Z a 



(84) 
(85) 
Az 2 a 2 . 

(86) 



therefore, s can be obtained by solving a simple nonlinear 
equation numerically 

The algorithm presented can be used to solve this 
one-dimensional problem. The planes are placed in the 
periodic box of size l\ x l 2 x Z 3 , at x\ = ±a. The 
amount of fixed charges, Ze, is distributed homoge- 
neously in the planes, so that the surface charge den- 
sity is cr = Ze/(2Z2^3)- The corresponding counterions 
are first distributed homogeneously in the region between 
the planes. The simulation box in x\ direction is some- 
what larger than 2a, in order to be able to implement 



8 




0.01 



FIG. 4: Counterion concentration profile along the x-axis 
for various grid resolutions as indicated by the legend. The 
continuous line is the solution of the Id isotropic cell model. 
The dash-dotted line shows the concentration profile for a 2d 
charged plate with the same surface charge density. 



periodic boundary conditions in this direction. Beyond 
the charged planes, both the concentration and the elec- 
tric field vanish. Due to the periodic boundary condi- 
tions in X2 and X3 direction, the system is translationally 
invariant in these directions, and hence the solution is 
one-dimensional. 

The calculations were performed for h — h — h ~ 
32A , a = 15A , Z = 60, and different grid sizes. The 
agreement between the simulation and the analytic ex- 
pression is quite good, see Fig. [TJ and increases with the 
grid resolution in x direction. Increasing the grid reso- 
lution in the orthogonal directions has no effect, as ex- 
pected (see Fig. (TJ lower part). Figure [2] shows the rela- 
tive error for different resolutions, denned by 



relative error = 



c a (x) 



Ca(x) 



(87) 



where c a (x) denotes the counterion concentration from 
the analytic solution and c num (x) is the numerical result. 



B. Colloidal Particle in a Box 

In this case a colloidal particle, modeled by a sphere 
of radius r, is placed at the center of the periodic box of 
size l\ x I2 x I3. The sphere carries a total charge of Ze, 
uniformly distributed over its surface. In the simulation, 
the fixed charges must be interpolated onto the nodes of 
the grid. The simplest way is to generate M»l random 
points distributed uniformly over the sphere surface. For 
each point, a (volume) charge density of Ze/(MAV) is 
added to the closest grid node. These nodes are marked 
as fixed. 

The calculations were performed for a cubic box of 
sides ii = I2 = fa = 30Aq. Grids with cubic symmetry 



T3 





£0.005- 



o 




□ x-axis 
o diagonal 




20 



FIG. 5: Electric field E along the a;-axis and along the (111) 
diagonal for a grid resolution of 256 3 . 



and various resolutions were used. A colloidal sphere 
of radius r = 3Ao and valence Z = 60 was placed at 
the center of the box, and the counterions were initially 
distributed uniformly over the outer space. 

Two independent versions of the algorithm were imple- 
mented, one in C++ |20j and another one in C. Runs for 
grids of different sizes were done on an Intel Core 2 Duo 
E6600 FSB 1066 2x2.4 Ghz, with 4GB RAM. The algo- 
rithm was run until the change in the functional reached 
a value smaller than 10~ 8 . Performance results are sum- 
marized in Table[U 

The pictures in Fig. [3] show a two-dimensional cut of 
the counterion concentration in a plane perpendicular to 
the z-axis, while in Fig. [4] one-dimensional cuts (along 
the ^-direction of the simulation box) are shown. The 
data show that near the surface of the charged colloid the 
discretization effects due to the cubic grid are quite large; 
nevertheless the data for the 256 3 grid seem reasonably 
well converged to the continuum limit. The concentra- 
tion profiles in (100) direction and in (111) direction are 
very close (i. e. within the resolution of the plot the 
curves coincide — these data have however not been in- 
cluded for the sake of clarity of the plot). Therefore one 
should expect that the solution is essentially identical 
to one obtained with strict spherical symmetry. This is 
indeed the case, as a comparison of the concentration 
profile with the corresponding solution of the spherically 
symmetric Poisson-Boltzmann cell model shows (see Fig. 
gj. 

In the latter, the cubic simulation cell is replaced by 
a spherical cell of the same volume, and the Poisson- 
Boltzmann equation is solved for the radial coordinate. 
In our calculations, this was done by transforming Eq. [5] 
to spherical coordinates, and then to a set of two coupled 
first-order differential equations (one for the potential, 
one for the field) . This set was solved by a simple integra- 
tor analogous to the velocity Verlet scheme known from 
Molecular Dynamics [2l[, using a step size of Ar = 10~ 4 
(in reduced units), and integrating from the colloid ra- 
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Code A (C++) 


Code A (C++) 


Code B (C) 


Code B (C) 


Grid Resolution 


Time (s) 


Memory (%) 


Time (s) 


Memory (%) 


32 x 32 x 32 


30 


0.2 


30 


0.1 


64 x 64 x 64 


836 


0.5 


844 


0.5 


128 x 128 x 128 


20701 


2.9 


20574 


3.6 


256 x 256 x 256 


381371 


22.5 


389030 


24.4 




Iterations 


Functional 


Iterations 


Functional 


32 x 32 x 32 


642 


1114.83 


634 


1113.88 


64 x 64 x 64 


2265 


1024.14 


2226 


1024.19 


128 x 128 x 128 


7206 


957.63 


7038 


957.63 


256 x 256 x 256 


20323 


923.85 


19711 


924.70 



TABLE I: Performance data for our two basic implementations. 



dius outwards. The constant A appearing in Eq. [5] was 
determined self-consistently by a shooting procedure, us- 
ing the requirement that the electric field must vanish at 
the outer radius, as a result of Gauss' law and the over- 
all charge neutrality. The thus-determined profile agrees 
quite well with the one obtained from our algorithm for 
the finest resolution. 

Nevertheless, the solution for the cubic geometry does 
exhibit some anisotropy that, by construction, is absent 
in the cell model. This is essentially invisible in the con- 
centration profile, but clearly observable in the electric 
field profile, as shown in Fig. [51 where the decays in 
(100) and (111) direction are compared. 

In the immediate vicinity of the colloidal surface, one 
may view the geometry as effectively planar. For a pla- 
nar surface, the solution is characterized by the so-called 
Gouy-Chapman length [22j. In our reduced units, this 
length has the value 0.044, which should be compared to 
the colloid radius (0.5728), and the lattice spacing (0.045 
for the 128 3 grid). The planar solution is shown in Fig. 
|4] (dash-dotted line); it also agrees reasonably well with 
the profiles from our algorithm. Altogether, the data 
indicate that a lattice spacing of roughly half a Gouy- 
Chapman length is small enough to yield a reasonably 
well converged solution. 



IV. SPEEDUPS 

As one sees from Tab. [Q the number of necessary iter- 
ations and the amount of CPU time are quite large. We 
have therefore looked for strategies to speed up the pro- 
cedure without sacrificing the basic formulation that pro- 
vides intrinsic stability. One possibility is to remove the 
rotational component of the electric field not by means 
of plaquette moves but rather by solving the Poisson 
equation. Within the chosen discretization scheme, the 
electrostatic potential needs to be an object associated 
with the sites, such that the electric field on a link is 
obtained by simply taking the potential difference be- 
tween the adjacent sites. This leads to the simplest pos- 



sible finite-difference scheme for the Poisson equation, 
which can be solved efficiently and in a stable way by 
using a Fast Fourier Transform (FFT) and the appropri- 
ate lattice Green's function [23(. This has the advantage 
that one single lattice sweep not only reduces the rota- 
tional component (as is the case for the plaquette moves) 
but rather eliminates it completely. Therefore the FFT 
promises to increase the convergence speed. An easy im- 
plementation is possible u sing the well-known and effi- 
cient fftw3 library routine 24j. It should be noted that 
the link moves, which update the concentrations and the 
fields simultaneously, remain unchanged, such that the 
procedure still stays strictly on the constraint surface. 

In a first implementation, we eliminated all plaquette 
moves and replaced them with FFT sweeps done dur- 
ing initialization as well as subsequently after every 25th 
link sweep. As seen from Tab. [Til this improves the effi- 
ciency roughly by a factor of 1.1 ... 2. These results were 
obtained on the same computer as those of Tab. Q] 

Furthermore, we can tackle the slowdown that comes 
from the fact that the ions have to be moved by site-by- 
site hops throughout the system ( "hydrodynamic slowing 
down"). To this end, we first run the calculation on a 
rather coarse grid (in practice, we started with 8x8x8), 
such that most of the necessary "mass transport" is al- 
ready done in that preliminary run. Starting from there, 
we go to a finer grid (in practice, we reduced the lat- 
tice spacing in all three directions by a factor of two) 
and linearly interpolate the output of the previous run 
onto that grid. Then the free energy is relaxed again; the 
output of that run is interpolated onto a yet finer grid, 
and so on. Obviously, this can be done rather easily by 
a straightforward recursion, until the desired grid reso- 
lution is reached. The runs before the finest resolution 
may then be viewed as a "pre-conditioner" . This opti- 
mization yields another speedup by roughly 25%, as seen 
from Tab. El 

These are obviously two rather simple optimizations, 
which do not interfere with the basic data structure of 
the simple Cartesian grid. Further optimizations, which 
are however much more complicated, are possible by (i) 
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Pure FFT 


Pure FFT 


FFT + pro-conditioner 


FFT + pre-conditioner 


Grid Resolution 


Time (s) 


Number of iterations 


Time (s) 


Number of iterations 


32 x 32 x 32 


14 


297 


11 


242 


64 x 64 x 64 


409 


1150 


331 


888 


128 x 128 x 128 


12180 


4427 


9744 


3421 


256 x 256 x 256 


354783 


16766 


273538 


12430 



TABLE II: Performance data for our two speeded-up implementations. Note that the data for the runs with pre-conditioner 
mean (i) CPU time for the overall procedure, and (ii) number of iterations in the final run with the finest resolution. 



adaptive mesh refinement (i. e. a fine resolution is only 
used in those regions where the fields vary strongly) , and 
(ii) using finite-element-type unstructured grids. Such 
more advanced approaches would be based upon con- 
structing the dual (Voronoi) lattice in order to define and 
calculate the fluxes. While this is expected to yield fur- 
ther substantial speedups, this was not attempted here, 
and is rather mentioned as a suggestion for the larger 
community. 

V. CONCLUDING REMARKS 

Variational techniques were used to develop an algo- 
rithm for the Poisson-Boltzmann equation. The required 
amount of memory scales linearly with the system size. 
In our first simple implementation, all moves are local. 
This leads to fast memory access and fast calculations 
at each move. Furthermore, the algorithm can be easily 
parallelized. On the other hand, the charges need time 
to move between distant nodes. In fact, while the time 
spent on each move is essentially constant, the total num- 
ber of sweeps required depends highly on how large a grid 
is swept. Considerable speedups were possible by using 
FFTs and a hierarchical pre-conditioner, but the basic 
"hydrodynamic slowing down" remains still present. Fur- 
ther speedups are likely to be possible by adaptive mesh 
refinement and by using unstructured grids. 

Existing algorithms based on standard discretizations 
of the differential operators are probably significantly 



faster than even the fastest of our current implementa- 
tions. Nevertheless, the approach that we have outlined 
in the present paper has the inherent advantage that 
its mathematical formulation provides intrinsic stability. 
This is mainly due to the fact that the free energy func- 
tional is formulated in terms of the electric field instead 
of the electrostatic potential, which provides a search for 
a true minimum in function space, rather than for a sad- 
dle point. As a result, the algorithm is very robust, in 
the sense that it will always converge to the solution, and 
at every iteration it approaches the solution more closely. 
Furthermore, the essential conservation laws that are at 
the heart of electrostatics — conservation of mass and 
conservation of electric flux — are built into the formu- 
lation with machine accuracy. We believe that these are 
fundamental advantages directly related to the under- 
lying physics of the problem, and we consider them as 
much more important than implementation details. The 
Maggs formulation provides a new way of thinking about 
electrostatics, and we hope that this, in combination with 
existing numerical "technology" , will in the future bring 
about very useful algorithms. 
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